Project: 1st Order Spatial Point Patterns Analysis Methods

Published

March 23, 2023

Modified

April 7, 2023

pacman::p_load(readxl, httr, jsonlite, maptools, sf, sfdep, raster, spatstat, spNetwork, rgdal, sp, tmap, tidyverse)

Geospatial Data Wrangling

For 1st order SPPA: airbnb_sf, mpsz_sf, sg_sf For 2nd order SPPA: same as above For NetSPPA: roadnetwork_sf For K function, K cross function, LCLQ: sevenele_sf, busstop_sf, hotel_sf, mrt_sf, mall_sf, attr_sf, unis_sf

read the 6 shp, kml, geojson files

busstop_sf <- st_read(dsn = "data/geospatial/busstop-022023-shp", layer="BusStop")
Reading layer `BusStop' from data source 
  `C:\valtyl\IS415-GAA\Project\data\geospatial\busstop-022023-shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 5159 features and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 3970.122 ymin: 26482.1 xmax: 48284.56 ymax: 52983.82
Projected CRS: SVY21
mpsz_sf <- st_read(dsn = "data/geospatial/MPSZ-2019-shp", layer="MPSZ-2019")
Reading layer `MPSZ-2019' from data source 
  `C:\valtyl\IS415-GAA\Project\data\geospatial\MPSZ-2019-shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 332 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 103.6057 ymin: 1.158699 xmax: 104.0885 ymax: 1.470775
Geodetic CRS:  WGS 84
network_sf <- st_read(dsn = "data/geospatial/roadnetwork-022023-shp", layer="RoadSectionLine")
Reading layer `RoadSectionLine' from data source 
  `C:\valtyl\IS415-GAA\Project\data\geospatial\roadnetwork-022023-shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 15300 features and 2 fields
Geometry type: MULTILINESTRING
Dimension:     XY
Bounding box:  xmin: 3248.394 ymin: 22755.1 xmax: 50043.6 ymax: 50135.95
Projected CRS: SVY21
attr_sf <- st_read(dsn = "data/geospatial/tourism-shp", layer="TOURISM")
Reading layer `TOURISM' from data source 
  `C:\valtyl\IS415-GAA\Project\data\geospatial\tourism-shp' using driver `ESRI Shapefile'
Simple feature collection with 109 features and 16 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 11380.23 ymin: 22858.67 xmax: 43659.54 ymax: 47627.69
Projected CRS: SVY21
mrt_sf <- st_read(dsn = "data/geospatial/mrt-112022-shp", layer="Train_Station_Exit_Layer") 
Reading layer `Train_Station_Exit_Layer' from data source 
  `C:\valtyl\IS415-GAA\Project\data\geospatial\mrt-112022-shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 562 features and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 6134.086 ymin: 27499.7 xmax: 45356.36 ymax: 47865.92
Projected CRS: SVY21
hotel_sf <- st_read("data/geospatial/hotels-2021-kml/hotel-locations.kml")
Reading layer `HOTELS' from data source 
  `C:\valtyl\IS415-GAA\Project\data\geospatial\hotels-2021-kml\hotel-locations.kml' 
  using driver `KML'
Simple feature collection with 422 features and 2 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 103.6351 ymin: 1.245797 xmax: 103.9891 ymax: 1.419278
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
sg_sf <- st_read(dsn = "data/geospatial/costal-shp", layer="CostalOutline")
Reading layer `CostalOutline' from data source 
  `C:\valtyl\IS415-GAA\Project\data\geospatial\costal-shp' using driver `ESRI Shapefile'
Simple feature collection with 60 features and 4 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 2663.926 ymin: 16357.98 xmax: 56047.79 ymax: 50244.03
Projected CRS: SVY21

read mall csv and convert to sf

mall_csv <- read_csv("data/geospatial/shoppingmalls-2019-csv/mall_coordinates_updated.csv")
mall_sf <- st_as_sf(mall_csv, coords = c("longitude", "latitude"), crs=4326)
mall_sf
Simple feature collection with 184 features and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 103.6784 ymin: 1.263797 xmax: 103.9897 ymax: 1.448227
Geodetic CRS:  WGS 84
# A tibble: 184 × 3
    ...1 name                               geometry
 * <dbl> <chr>                           <POINT [°]>
 1     0 100 AM                  (103.8435 1.274588)
 2     1 112 KATONG              (103.9051 1.305087)
 3     2 313@SOMERSET            (103.8377 1.301385)
 4     3 321 CLEMENTI             (103.765 1.312025)
 5     4 600 @ TOA PAYOH          (103.851 1.334042)
 6     5 888 PLAZA               (103.7953 1.437131)
 7     6 AMK HUB                 (103.8485 1.369223)
 8     7 ADMIRALTY PLACE         (103.8018 1.439881)
 9     8 ALEXANDRA RETAIL CENTRE (103.8014 1.273843)
10     9 ANCHORPOINT             (103.8056 1.288935)
# … with 174 more rows

read airbnb csv and convert to sf

airbnb_csv <- read_csv("data/geospatial/airbnb-2022-csv/listings.csv")
airbnb_sf <- st_as_sf(airbnb_csv, coords = c("longitude", "latitude"), crs=4326)
get_coords <- function(add_list){
  
  # Create a data frame to store all retrieved coordinates
  postal_coords <- data.frame()
  
  # loop to go through each address in the list  
  for (i in add_list){
    #print(i)
    
    # response from OneMap API
    r <- GET('https://developers.onemap.sg/commonapi/search?',
           query=list(searchVal=i,
                     returnGeom='Y',
                     getAddrDetails='Y'))
    data <- fromJSON(rawToChar(r$content))
    found <- data$found
    res <- data$results
    
    # Create a new data frame for each address
    new_row <- data.frame()
    
    # If single result, append 
    if (found == 1){
      postal <- res$POSTAL 
      lat <- res$LATITUDE
      lng <- res$LONGITUDE
      new_row <- data.frame(address= i, postal = postal, latitude = lat, longitude = lng)
    }
    
    # If multiple results, drop NIL and append top 1
    else if (found > 1){
      # Remove those with NIL as postal
      res_sub <- res[res$POSTAL != "NIL", ]
      
      # Set as NA first if no Postal
      if (nrow(res_sub) == 0) {
          new_row <- data.frame(address= i, postal = NA, latitude = NA, longitude = NA)
      }
      
      else{
        top1 <- head(res_sub, n = 1)
        postal <- top1$POSTAL 
        lat <- top1$LATITUDE
        lng <- top1$LONGITUDE
        new_row <- data.frame(address= i, postal = postal, latitude = lat, longitude = lng)
      }
    }

    else {
      new_row <- data.frame(address= i, postal = NA, latitude = NA, longitude = NA)
    }
    
    # Add the row to the main dataframe
    postal_coords <- rbind(postal_coords, new_row)
  }
  return(postal_coords)
}

read 7-11 xlsx and get coords

Show the code!
sevenele_xlsx <- read_xlsx("data/geospatial/711stores-032023-xlsx/7-11_convenience_stores.xlsx")

sevenele_list <- sort(unique(sevenele_xlsx$Postal))
sevenele_coords <- get_coords(sevenele_list)
sevenele_coords[(is.na(sevenele_coords$postal) | is.na(sevenele_coords$latitude) | is.na(sevenele_coords$longitude) | sevenele_coords$postal=="NIL"), ]

fix 711 coords

Show the code!
sevenele_coords[64,]$postal <- "098975"
sevenele_coords[64,]$latitude <- "1.255169"
sevenele_coords[64,]$longitude <- "103.818508"

sevenele_coords[455,]$postal <- "828694"
sevenele_coords[455,]$latitude <- "1.420668"
sevenele_coords[455,]$longitude <- "103.912257"

check if fixed

Show the code!
sevenele_coords[(is.na(sevenele_coords$postal) | is.na(sevenele_coords$latitude) | is.na(sevenele_coords$longitude) | sevenele_coords$postal=="NIL"), ]

get 711 rds

Show the code!
sevenele_df <- left_join(sevenele_xlsx, sevenele_coords, by=c('Postal' = 'address'))

sevenele_rds <- write_rds(sevenele_df, "data/geospatial/711stores-032023-xlsx/sevenele_rds.rds")

get 711 sf

seveneles <- read_rds("data/geospatial/711stores-032023-xlsx/sevenele_rds.rds")
sevenele_sf <- st_as_sf(seveneles,
                    coords = c("longitude", 
                               "latitude"),
                    crs=4326) %>%
  st_transform(crs = 3414)

read universities and colleges xlsx and get coords

Show the code!
unis_xlsx <- read_xlsx("data/geospatial/universitiescolleges-2023-xlsx/universities_and_colleges.xlsx")

unis_list <- sort(unique(unis_xlsx$name))
unis_coords <- get_coords(unis_list)
unis_coords[(is.na(unis_coords$postal) | is.na(unis_coords$latitude) | is.na(unis_coords$longitude) | unis_coords$postal=="NIL"), ]

fix uni coords

Show the code!
unis_coords[1,]$postal <- "179104"
unis_coords[1,]$latitude <- "1.29006"
unis_coords[1,]$longitude <- "103.8497"

unis_coords[2,]$postal <- "229469"
unis_coords[2,]$latitude <- "1.31648392303"
unis_coords[2,]$longitude <- "103.873964955"

unis_coords[4,]$postal <- "408601"
unis_coords[4,]$latitude <- "1.3196562"
unis_coords[4,]$longitude <- "103.8923565"

unis_coords[5,]$postal <- "059405"
unis_coords[5,]$latitude <- "1.28799"
unis_coords[5,]$longitude <- "103.84692"

unis_coords[6,]$postal <- "188946"
unis_coords[6,]$latitude <- "1.3462717"
unis_coords[6,]$longitude <- "103.95302009"

unis_coords[7,]$postal <- "329162"
unis_coords[7,]$latitude <- "1.3270276324741"
unis_coords[7,]$longitude <- "103.85336279869"

unis_coords[10,]$postal <- "188647"
unis_coords[10,]$latitude <- "1.2990616"
unis_coords[10,]$longitude <- "103.8494881"

unis_coords[11,]$postal <- "088383"
unis_coords[11,]$latitude <- "1.28013"
unis_coords[11,]$longitude <- "103.84193"

unis_coords[12,]$postal <- "049319"
unis_coords[12,]$latitude <- "1.28469"
unis_coords[12,]$longitude <- "103.85269"

unis_coords[14,]$postal <- "138676"
unis_coords[14,]$latitude <- "1.299122"
unis_coords[14,]$longitude <- "103.788025"

unis_coords[16,]$postal <- "228095"
unis_coords[16,]$latitude <- "1.30209070954"
unis_coords[16,]$longitude <- "103.849560613"

unis_coords[17,]$postal <- "188306"
unis_coords[17,]$latitude <- "1.3001799"
unis_coords[17,]$longitude <- "103.84924"

unis_coords[19,]$postal <- "069542"
unis_coords[19,]$latitude <- "1.2750794"
unis_coords[19,]$longitude <- "103.8462627"

unis_coords[20,]$postal <- "148951"
unis_coords[20,]$latitude <- "1.2970726"
unis_coords[20,]$longitude <- "103.8011673"

unis_coords[22,]$postal <- "189655"
unis_coords[22,]$latitude <- "1.30020700501"
unis_coords[22,]$longitude <- "103.851710836"

unis_coords[26,]$postal <- "248922"
unis_coords[26,]$latitude <- "1.31532"
unis_coords[26,]$longitude <- "103.77994"

unis_coords[30,]$postal <- "599491"
unis_coords[30,]$latitude <- "1.326"
unis_coords[30,]$longitude <- "103.79999"

unis_coords[37,]$postal <- "550266"
unis_coords[37,]$latitude <- "1.35324"
unis_coords[37,]$longitude <- "103.87145"

unis_coords[38,]$postal <- "139660"
unis_coords[38,]$latitude <- "1.3079"
unis_coords[38,]$longitude <- "103.777"

check fix

Show the code!
unis_coords[(is.na(unis_coords$postal) | is.na(unis_coords$latitude) | is.na(unis_coords$longitude) | unis_coords$postal=="NIL"), ]

get uni rds

Show the code!
unis_df <- left_join(unis_xlsx, unis_coords, by=c('name' = 'address'))
unis_rds <- write_rds(unis_df, "data/geospatial/universitiescolleges-2023-xlsx/unis_rds.rds")

get unis sf

unis <- read_rds("data/geospatial/universitiescolleges-2023-xlsx/unis_rds.rds")
unis_sf <- st_as_sf(unis,
                    coords = c("longitude", 
                               "latitude"),
                    crs=4326) %>%
  st_transform(crs = 3414)

remove unnecessary columns (except mpsz_sf, network_sf)

attr_sf <- attr_sf %>% select(c(5))
busstop_sf <- busstop_sf %>% select(c(1))
hotel_sf <- hotel_sf %>% select(c(1))
mall_sf <- mall_sf %>% select(c(2))
mrt_sf <- mrt_sf %>% select(c(1))
sevenele_sf <- sevenele_sf %>% select(c(1))
unis_sf <- unis_sf %>% select(c(1))
airbnb_sf <- airbnb_sf %>% select(c(2,6,7,8))

remove z dimensions

hotel_sf <- st_zm(hotel_sf)

check for invalid geometries

length(which(st_is_valid(attr_sf) == FALSE))
[1] 0
length(which(st_is_valid(busstop_sf) == FALSE))
[1] 0
length(which(st_is_valid(hotel_sf) == FALSE))
[1] 0
length(which(st_is_valid(mall_sf) == FALSE))
[1] 0
length(which(st_is_valid(mrt_sf) == FALSE))
[1] 0
length(which(st_is_valid(sevenele_sf) == FALSE))
[1] 0
length(which(st_is_valid(unis_sf) == FALSE))
[1] 0
length(which(st_is_valid(mpsz_sf) == FALSE))
[1] 6
length(which(st_is_valid(network_sf) == FALSE))
[1] 0
length(which(st_is_valid(airbnb_sf) == FALSE))
[1] 0
length(which(st_is_valid(sg_sf) == FALSE))
[1] 1
hotel_sf <- st_zm(hotel_sf)

fix mpsz invalid geometry and check again

mpsz_sf <- st_make_valid(mpsz_sf)
length(which(st_is_valid(mpsz_sf) == FALSE))
[1] 0

fix sg invalid geometry and check again

sg_sf <- st_make_valid(sg_sf)
length(which(st_is_valid(sg_sf) == FALSE))
[1] 0

dealing w missing values

attr_sf[rowSums(is.na(attr_sf))!=0,]
Simple feature collection with 0 features and 1 field
Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
Projected CRS: SVY21
[1] PAGETITLE geometry 
<0 rows> (or 0-length row.names)
busstop_sf[rowSums(is.na(busstop_sf))!=0,]
Simple feature collection with 0 features and 1 field
Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
Projected CRS: SVY21
[1] BUS_STOP_N geometry  
<0 rows> (or 0-length row.names)
hotel_sf[rowSums(is.na(hotel_sf))!=0,]
Simple feature collection with 0 features and 1 field
Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
Geodetic CRS:  WGS 84
[1] Name     geometry
<0 rows> (or 0-length row.names)
mall_sf[rowSums(is.na(mall_sf))!=0,]
Simple feature collection with 0 features and 1 field
Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
Geodetic CRS:  WGS 84
# A tibble: 0 × 2
# … with 2 variables: name <chr>, geometry <GEOMETRY [°]>
mrt_sf[rowSums(is.na(mrt_sf))!=0,]
Simple feature collection with 0 features and 1 field
Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
Projected CRS: SVY21
[1] stn_name geometry
<0 rows> (or 0-length row.names)
sevenele_sf[rowSums(is.na(sevenele_sf))!=0,]
Simple feature collection with 0 features and 1 field
Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
Projected CRS: SVY21 / Singapore TM
# A tibble: 0 × 2
# … with 2 variables: Name <chr>, geometry <GEOMETRY [m]>
unis_sf[rowSums(is.na(unis_sf))!=0,]
Simple feature collection with 0 features and 1 field
Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
Projected CRS: SVY21 / Singapore TM
# A tibble: 0 × 2
# … with 2 variables: name <chr>, geometry <GEOMETRY [m]>
mpsz_sf[rowSums(is.na(mpsz_sf))!=0,]
Simple feature collection with 0 features and 6 fields
Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
Geodetic CRS:  WGS 84
[1] SUBZONE_N  SUBZONE_C  PLN_AREA_N PLN_AREA_C REGION_N   REGION_C   geometry  
<0 rows> (or 0-length row.names)
network_sf[rowSums(is.na(network_sf))!=0,]
Simple feature collection with 15300 features and 2 fields
Geometry type: MULTILINESTRING
Dimension:     XY
Bounding box:  xmin: 3248.394 ymin: 22755.1 xmax: 50043.6 ymax: 50135.95
Projected CRS: SVY21
First 10 features:
   RD_CD            RD_CD_DESC                       geometry
1   <NA>       SEA BREEZE ROAD MULTILINESTRING ((41260.55 ...
2   <NA>       HILLVIEW AVENUE MULTILINESTRING ((20404.4 3...
3   <NA>     SERANGOON CENTRAL MULTILINESTRING ((32508.6 3...
4   <NA> CHANGI SOUTH STREET 3 MULTILINESTRING ((42949.75 ...
5   <NA>           PANDAN ROAD MULTILINESTRING ((18683.23 ...
6   <NA>      WEST COAST GROVE MULTILINESTRING ((19396.72 ...
7   <NA>          NEYTHAL ROAD MULTILINESTRING ((12917.95 ...
8   <NA>         WILLOW AVENUE MULTILINESTRING ((32638.54 ...
9   <NA> CHOA CHU KANG NORTH 5 MULTILINESTRING ((18182.78 ...
10  <NA>       CASSIA CRESCENT MULTILINESTRING ((33711.34 ...
airbnb_sf[rowSums(is.na(airbnb_sf))!=0,]
Simple feature collection with 0 features and 4 fields
Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
Geodetic CRS:  WGS 84
# A tibble: 0 × 5
# … with 5 variables: name <chr>, neighbourhood <chr>, room_type <chr>,
#   price <dbl>, geometry <GEOMETRY [°]>
sg_sf[rowSums(is.na(sg_sf))!=0,]
Simple feature collection with 0 features and 4 fields
Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
Projected CRS: SVY21
[1] GDO_GID    MSLINK     MAPID      COSTAL_NAM geometry  
<0 rows> (or 0-length row.names)

verify coordinate system

st_crs(attr_sf)
Coordinate Reference System:
  User input: SVY21 
  wkt:
PROJCRS["SVY21",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["unnamed",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]
st_crs(busstop_sf)
Coordinate Reference System:
  User input: SVY21 
  wkt:
PROJCRS["SVY21",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["unnamed",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]
st_crs(hotel_sf)
Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]
st_crs(mall_sf)
Coordinate Reference System:
  User input: EPSG:4326 
  wkt:
GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]
st_crs(mrt_sf)
Coordinate Reference System:
  User input: SVY21 
  wkt:
PROJCRS["SVY21",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["unnamed",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]
st_crs(sevenele_sf)
Coordinate Reference System:
  User input: EPSG:3414 
  wkt:
PROJCRS["SVY21 / Singapore TM",
    BASEGEOGCRS["SVY21",
        DATUM["SVY21",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4757]],
    CONVERSION["Singapore Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping."],
        AREA["Singapore - onshore and offshore."],
        BBOX[1.13,103.59,1.47,104.07]],
    ID["EPSG",3414]]
st_crs(unis_sf)
Coordinate Reference System:
  User input: EPSG:3414 
  wkt:
PROJCRS["SVY21 / Singapore TM",
    BASEGEOGCRS["SVY21",
        DATUM["SVY21",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4757]],
    CONVERSION["Singapore Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping."],
        AREA["Singapore - onshore and offshore."],
        BBOX[1.13,103.59,1.47,104.07]],
    ID["EPSG",3414]]
st_crs(mpsz_sf)
Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]
st_crs(network_sf)
Coordinate Reference System:
  User input: SVY21 
  wkt:
PROJCRS["SVY21",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["unnamed",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]
st_crs(airbnb_sf)
Coordinate Reference System:
  User input: EPSG:4326 
  wkt:
GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]
st_crs(sg_sf)
Coordinate Reference System:
  User input: SVY21 
  wkt:
PROJCRS["SVY21",
    BASEGEOGCRS["SVY21[WGS84]",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["unnamed",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]

transform coordinate system

# with st_set_crs(), we can assign the appropriate ESPG Code
attr_sf <- st_set_crs(attr_sf, 3414)
busstop_sf <- st_set_crs(busstop_sf, 3414)
mrt_sf <- st_set_crs(mrt_sf, 3414)
network_sf <- st_set_crs(network_sf, 3414)
sg_sf <- st_set_crs(sg_sf, 3414)

# with st_transform(), we can change from one CRS to another
mpsz_sf <- st_transform(mpsz_sf, crs=3414)
hotel_sf <- st_transform(hotel_sf, crs=3414)
mall_sf <- st_transform(mall_sf, crs=3414)
airbnb_sf <- st_transform(airbnb_sf, crs=3414)

# sevenele_sf and unis_sf are in the correct CRS and ESPG code

Preparing for 1st Order SPPA

convert sf to sp’s spatial class

airbnb <- as_Spatial(airbnb_sf)
mpsz <- as_Spatial(mpsz_sf)
sg <- as_Spatial(sg_sf)

convert sp’s spatial into generic sp format

airbnb_sp <- as(airbnb, "SpatialPoints")
sg_sp <- as(sg, "SpatialPolygons")
airbnb
class       : SpatialPointsDataFrame 
features    : 3037 
extent      : 7406.989, 44064, 25352.92, 48321.84  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
variables   : 4
names       :                                               name, neighbourhood,       room_type,  price 
min values  :      -Master bedroom with attached bathroom at CBD,    Ang Mo Kio, Entire home/apt,      0 
max values  : 阿裕尼公寓普通房急需女搭房一名联系人:郑小姐 电话:,        Yishun,     Shared room, 135140 

convert generic sp format into spatstat’s ppp format

airbnb_ppp <- as(airbnb_sp, "ppp")
airbnb_ppp
Planar point pattern: 3037 points
window: rectangle = [7406.99, 44064] x [25352.92, 48321.84] units

check number of duplicate point events

sum(multiplicity(airbnb_ppp) > 1)
[1] 211

handle duplicate points:

airbnb_ppp_jit <- rjitter(airbnb_ppp, 
                             retry=TRUE, 
                             nsim=1, 
                             drop=TRUE)

check duplicated points:

any(duplicated(airbnb_ppp_jit))
[1] FALSE

create owin

# sg_owin <- as(sg_sp, "owin")
central = mpsz_sf[mpsz_sf$REGION_N=='CENTRAL REGION',]
central_owin <- central %>%
  as ('Spatial') %>%
  as ('SpatialPolygons') %>%
  as ('owin')
plot(central_owin)

combine airbnb point events obj and owin obj

airbnbSG_ppp = airbnb_ppp_jit[central_owin]

summary

summary(airbnbSG_ppp)
Planar point pattern:  2476 points
Average intensity 1.815334e-05 points per square unit

Coordinates are given to 3 decimal places
i.e. rounded to the nearest multiple of 0.001 units

Window: polygonal boundary
164 separate polygons (14 holes)
                  vertices         area relative.area
polygon 1              299  1.84404e+06      1.35e-02
polygon 2              164  3.92563e+05      2.88e-03
polygon 3              238  5.06591e+05      3.71e-03
polygon 4              132  3.88733e+05      2.85e-03
polygon 5              255  1.59034e+06      1.17e-02
polygon 6             1020  1.27781e+06      9.37e-03
polygon 7 (hole)         4 -1.23063e-04     -9.02e-13
polygon 8 (hole)         3 -1.06327e-01     -7.80e-10
polygon 9 (hole)         5 -5.18342e-01     -3.80e-09
polygon 10 (hole)        3 -2.53510e-02     -1.86e-10
polygon 11 (hole)        4 -7.56223e-02     -5.54e-10
polygon 12 (hole)        4 -7.19058e-01     -5.27e-09
polygon 13             211  4.70522e+05      3.45e-03
polygon 14             155  2.67503e+05      1.96e-03
polygon 15             129  9.53762e+04      6.99e-04
polygon 16              94  5.96175e+04      4.37e-04
polygon 17              59  3.43163e+04      2.52e-04
polygon 18              10  4.90997e+02      3.60e-06
polygon 19               6  4.50376e+02      3.30e-06
polygon 20               4  2.69408e+02      1.98e-06
polygon 21            1424  4.87153e+06      3.57e-02
polygon 22 (hole)        4 -2.46124e-04     -1.80e-12
polygon 23 (hole)        3 -3.69186e-03     -2.71e-11
polygon 24 (hole)        3 -2.73199e-02     -2.00e-10
polygon 25 (hole)       11 -8.36614e+01     -6.13e-07
polygon 26              75  1.73525e+04      1.27e-04
polygon 27              40  1.38603e+04      1.02e-04
polygon 28              83  5.28926e+03      3.88e-05
polygon 29             137  3.22285e+03      2.36e-05
polygon 30             147  3.10396e+03      2.28e-05
polygon 31             106  3.04131e+03      2.23e-05
polygon 32              45  2.51228e+03      1.84e-05
polygon 33             438  3.45157e+06      2.53e-02
polygon 34              83  1.03237e+05      7.57e-04
polygon 35             102  1.12730e+06      8.27e-03
polygon 36             829  2.70003e+06      1.98e-02
polygon 37 (hole)       23 -1.99689e+01     -1.46e-07
polygon 38 (hole)       35 -1.38430e+02     -1.01e-06
polygon 39 (hole)       19 -4.37731e+00     -3.21e-08
polygon 40 (hole)      270 -1.21366e+03     -8.90e-06
polygon 41              88  5.33462e+05      3.91e-03
polygon 42              95  1.45518e+05      1.07e-03
polygon 43              55  6.35706e+05      4.66e-03
polygon 44              53  2.76828e+05      2.03e-03
polygon 45             114  6.36655e+04      4.67e-04
polygon 46              83  1.96620e+05      1.44e-03
polygon 47              33  3.65334e+05      2.68e-03
polygon 48             110  1.45503e+06      1.07e-02
polygon 49             135  8.53202e+05      6.26e-03
polygon 50             196  1.07072e+06      7.85e-03
polygon 51              47  5.33012e+05      3.91e-03
polygon 52              85  4.42296e+05      3.24e-03
polygon 53              38  4.11723e+05      3.02e-03
polygon 54             227  5.87221e+05      4.31e-03
polygon 55              35  3.94370e+04      2.89e-04
polygon 56              96  1.88768e+05      1.38e-03
polygon 57              59  1.33007e+05      9.75e-04
polygon 58              46  4.48128e+05      3.29e-03
polygon 59              31  5.21199e+05      3.82e-03
polygon 60              17  3.50788e+05      2.57e-03
polygon 61              54  2.61841e+05      1.92e-03
polygon 62             152  1.63038e+06      1.20e-02
polygon 63             180  5.61608e+05      4.12e-03
polygon 64              47  1.60809e+05      1.18e-03
polygon 65              49  5.95241e+05      4.36e-03
polygon 66              49  3.87611e+05      2.84e-03
polygon 67              59  1.03038e+06      7.55e-03
polygon 68              83  5.51731e+05      4.05e-03
polygon 69              69  2.90185e+05      2.13e-03
polygon 70             217  1.08479e+06      7.95e-03
polygon 71              41  6.28892e+05      4.61e-03
polygon 72             226  1.82685e+06      1.34e-02
polygon 73              55  2.93699e+05      2.15e-03
polygon 74             247  5.57280e+05      4.09e-03
polygon 75              48  5.56817e+04      4.08e-04
polygon 76              60  1.16331e+05      8.53e-04
polygon 77             343  2.00345e+06      1.47e-02
polygon 78             129  2.43459e+06      1.78e-02
polygon 79              58  3.10535e+05      2.28e-03
polygon 80             114  1.38034e+06      1.01e-02
polygon 81             134  1.95176e+06      1.43e-02
polygon 82             270  4.51538e+05      3.31e-03
polygon 83             123  1.72655e+05      1.27e-03
polygon 84             137  3.36218e+05      2.47e-03
polygon 85              62  7.42203e+05      5.44e-03
polygon 86             114  1.07899e+06      7.91e-03
polygon 87             319  4.60551e+05      3.38e-03
polygon 88             198  5.43484e+05      3.98e-03
polygon 89              51  2.78303e+05      2.04e-03
polygon 90             297  8.86957e+05      6.50e-03
polygon 91             182  2.18809e+05      1.60e-03
polygon 92             129  2.00787e+05      1.47e-03
polygon 93             169  7.10568e+05      5.21e-03
polygon 94              34  7.48683e+05      5.49e-03
polygon 95             173  3.68483e+05      2.70e-03
polygon 96             294  7.60622e+06      5.58e-02
polygon 97             238  2.21999e+05      1.63e-03
polygon 98             130  2.80176e+05      2.05e-03
polygon 99             140  2.14251e+05      1.57e-03
polygon 100             83  1.73123e+05      1.27e-03
polygon 101            192  5.91779e+05      4.34e-03
polygon 102            174  1.75348e+06      1.29e-02
polygon 103            192  3.40742e+05      2.50e-03
polygon 104            218  3.29438e+05      2.42e-03
polygon 105             87  1.70664e+05      1.25e-03
polygon 106            217  1.34616e+06      9.87e-03
polygon 107             27  1.71336e+05      1.26e-03
polygon 108             84  4.96255e+04      3.64e-04
polygon 109            198  1.93991e+05      1.42e-03
polygon 110             77  1.20171e+05      8.81e-04
polygon 111            273  6.14923e+05      4.51e-03
polygon 112            108  1.02656e+06      7.53e-03
polygon 113            154  1.67538e+05      1.23e-03
polygon 114             81  1.16002e+06      8.50e-03
polygon 115             86  9.63497e+05      7.06e-03
polygon 116            103  1.42508e+06      1.04e-02
polygon 117             74  1.22990e+06      9.02e-03
polygon 118             75  1.26341e+06      9.26e-03
polygon 119             50  3.69771e+05      2.71e-03
polygon 120             96  1.10727e+06      8.12e-03
polygon 121            126  3.38725e+06      2.48e-02
polygon 122             94  1.87800e+06      1.38e-02
polygon 123             40  8.67749e+05      6.36e-03
polygon 124             55  6.39144e+05      4.69e-03
polygon 125             54  4.11402e+05      3.02e-03
polygon 126             75  4.18656e+05      3.07e-03
polygon 127             74  2.18107e+05      1.60e-03
polygon 128            104  2.13519e+05      1.57e-03
polygon 129            119  4.85574e+05      3.56e-03
polygon 130            102  7.56923e+05      5.55e-03
polygon 131            117  3.51306e+05      2.58e-03
polygon 132             67  9.47451e+05      6.95e-03
polygon 133             93  1.05203e+06      7.71e-03
polygon 134             92  4.10940e+05      3.01e-03
polygon 135             72  8.39679e+05      6.16e-03
polygon 136            184  1.22706e+06      9.00e-03
polygon 137            120  5.58761e+05      4.10e-03
polygon 138            212  2.09608e+06      1.54e-02
polygon 139            119  2.15829e+06      1.58e-02
polygon 140            140  1.34746e+06      9.88e-03
polygon 141            381  2.03851e+06      1.49e-02
polygon 142            125  2.57843e+06      1.89e-02
polygon 143            112  7.35502e+05      5.39e-03
polygon 144            132  1.31911e+06      9.67e-03
polygon 145            122  1.37683e+06      1.01e-02
polygon 146            129  1.92662e+06      1.41e-02
polygon 147             74  1.79446e+06      1.32e-02
polygon 148            136  3.14493e+06      2.31e-02
polygon 149            195  2.63648e+06      1.93e-02
polygon 150             80  1.05717e+06      7.75e-03
polygon 151             69  4.39648e+05      3.22e-03
polygon 152             61  4.46241e+05      3.27e-03
polygon 153             72  5.72500e+05      4.20e-03
polygon 154             79  8.13382e+05      5.96e-03
polygon 155             97  1.03728e+06      7.61e-03
polygon 156            182  2.77329e+06      2.03e-02
polygon 157            131  3.46851e+06      2.54e-02
polygon 158             71  2.35998e+05      1.73e-03
polygon 159            244  4.70674e+05      3.45e-03
polygon 160             97  9.13532e+04      6.70e-04
polygon 161            128  1.71195e+06      1.26e-02
polygon 162            163  2.96147e+06      2.17e-02
polygon 163            268  3.84951e+06      2.82e-02
polygon 164            102  1.96415e+06      1.44e-02
enclosing rectangle: [18752.4, 37592.19] x [21678.34, 38891.36] units
                     (18840 x 17210 units)
Window area = 136394000 square units
Fraction of frame area: 0.421

1st Order SPPA

rescale to convert from m to km

airbnbSG_ppp.km <- rescale(airbnbSG_ppp, 1000, "km")

Kernel Density Estimation (automatic bandwidth)

kde_airbnbSG.bw <- density(airbnbSG_ppp.km,
                              sigma=bw.diggle,
                              edge=TRUE,
                            kernel="gaussian") 

other methods to do: bw.CvL(), bw.scott(), bw.ppl() other smoothing kernels: epanechnikov, quartic, disc

plot(kde_airbnbSG.bw)

KDE (fixed bandwidth)

kde_airbnbSG_600 <- density(airbnbSG_ppp.km, sigma=0.6, edge=TRUE, kernel="gaussian") #0.6 sigma = 600 m for bandwidth. 600m = 0.6km
plot(kde_airbnbSG_600)

sigma: range from 0.1 to 10 different kernels

KDE (adaptive bandwidth)

kde_airbnbSG_adaptive <- adaptive.density(airbnbSG_ppp.km, method="kernel")
plot(kde_airbnbSG_adaptive)

convert KDE output into grid (do for all)

gridded_kde_airbnbSG_bw <- as.SpatialGridDataFrame.im(kde_airbnbSG.bw)

# convert grid output into raster
kde_airbnbSG_bw_raster <- raster(gridded_kde_airbnbSG_bw)

# assign projection systems
projection(kde_airbnbSG_bw_raster) <- CRS("+init=EPSG:3414 +datum=WGS84 +units=km")
kde_airbnbSG_bw_raster
class      : RasterLayer 
dimensions : 128, 128, 16384  (nrow, ncol, ncell)
resolution : 0.1471859, 0.1344767  (x, y)
extent     : 18.7524, 37.59219, 21.67834, 38.89136  (xmin, xmax, ymin, ymax)
crs        : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +datum=WGS84 +units=km +no_defs 
source     : memory
names      : v 
values     : -3.377976e-13, 2071.377  (min, max)
kl = mpsz_sf[mpsz_sf$PLN_AREA_N == "KALLANG",]
dc = mpsz_sf[mpsz_sf$PLN_AREA_N == "DOWNTOWN CORE",]
ot = mpsz_sf[mpsz_sf$PLN_AREA_N == "OUTRAM",]
rc = mpsz_sf[mpsz_sf$PLN_AREA_N == "ROCHOR",]
mapex <- st_bbox(central)

final visualisation tmap

tmap_mode("plot")
tm_shape(kde_airbnbSG_bw_raster, bbox=mapex) + 
  tm_raster("v") +
  tm_layout(legend.position = c("right", "bottom"), frame = FALSE)

2nd Order SPPA G/F/K/L Function

G Function estimation

G_airbnb = Gest(airbnbSG_ppp, correction =c("rs", "Hanisch"))
plot(G_airbnb, xlim=c(0,500))

other correction: “best” Gest() - corrections: “none”, “rs”, “km”, “Hanisch”, “best”, “all” Fest() - corrections: “none”, “rs”, “km”, “cs”, “best”, “all” Kest() - corrections: “none”, “border”, “bord.modif”, “isotropic”, “Ripley”, “translate”, “translation”, “rigid”, “none”, “good”, “best”, “all” Lest() - same as Kest()

Complete Spatial Randomness Test

G_airbnb.csr <- envelope(airbnbSG_ppp, Gest, nsim = 999)
# G_airbnb.csr <- envelope(airbnbSG_ppp, Gest, correction="all", nsim = 999)

nsim: number btwn 0 to 1000

plot to show

plot(G_airbnb.csr)

K Function estimation

K_airbnb = Kest(airbnbSG_ppp, correction =c("best"))
plot(K_airbnb, xlim=c(0,500))

Complete Spatial Randomness Test

K_airbnb.csr <- envelope(airbnbSG_ppp, Kest, nsim = 999)
# G_airbnb.csr <- envelope(airbnbSG_ppp, Gest, correction="all", nsim = 999)

nsim: number btwn 0 to 1000

plot to show

plot(K_airbnb.csr)

visualisation

tmap_mode("plot")
tm_shape(mpsz_sf) + 
  tm_borders(alpha = 0.5) +
tm_shape(airbnb_sf) + 
  tm_dots(col = "red", size = 0.05) +
  tm_layout(main.title = "Airbnbs",
            main.title.position = 'center',
            main.title.size = 0.8,
            frame = TRUE)

occurrence <- airbnb_sf %>% count(airbnb_sf$neighbourhood)
occurrence <- occurrence[order(-occurrence$n),]
occurrence
Simple feature collection with 43 features and 2 fields
Geometry type: GEOMETRY
Dimension:     XY
Bounding box:  xmin: 7406.989 ymin: 25352.92 xmax: 44064 ymax: 48321.84
Projected CRS: SVY21 / Singapore TM
# A tibble: 43 × 3
   `airbnb_sf$neighbourhood`     n                                      geometry
   <chr>                     <int>                              <MULTIPOINT [m]>
 1 Kallang                     353 ((29508.14 32976.25), (29791.93 33174.18), (…
 2 Downtown Core               325 ((28945.03 28564.31), (28952.82 28557.68), (…
 3 Outram                      255 ((28290.63 29757.42), (28596.68 29346.08), (…
 4 Rochor                      214 ((29383.5 31768.77), (29521.51 31404.98), (2…
 5 Novena                      182 ((28059.14 33734.79), (28062.48 33523.59), (…
 6 Queenstown                  176 ((20006.09 30435.35), (20082.89 30453.04), (…
 7 Bukit Merah                 163 ((24405.42 28144.15), (24581.26 28085.54), (…
 8 River Valley                157 ((27291.23 31642.72), (27473.75 31495.65), (…
 9 Geylang                     145 ((32942.57 32745.19), (33084.45 32747.7), (3…
10 Bedok                       120 ((36040.87 33769.19), (36077.61 33347.9), (3…
# … with 33 more rows

kallang has the most airbnbs

K cross function

NetSPPA

kallang airbnbs:

kallang_airbnb_sf <- filter(airbnb_sf, neighbourhood == "Kallang")
kallang_airbnb_sf
Simple feature collection with 353 features and 4 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 29508.14 ymin: 31119.73 xmax: 33791.74 ymax: 34536.48
Projected CRS: SVY21 / Singapore TM
# A tibble: 353 × 5
   name                          neigh…¹ room_…² price            geometry
 * <chr>                         <chr>   <chr>   <dbl>         <POINT [m]>
 1 EnsuiteA w Balcony view/acce… Kallang Privat…    57 (31346.67 32390.22)
 2 Ensuite M1B within Heritage … Kallang Privat…    50    (30798 32708.67)
 3 EnsuiteK w Balcony view/acce… Kallang Privat…    57 (31239.83 32106.04)
 4 Heritage apartment: Master r… Kallang Privat…    69 (30505.31 32834.72)
 5 Lavender QueenBed *rmFr: 5m … Kallang Privat…    59 (31368.93 32169.07)
 6 Life Impact Coaching          Kallang Entire…   130 (32498.52 32598.12)
 7 1-Pax Small Private Room, Sh… Kallang Privat…    78 (31000.55 32666.66)
 8 2-pax Small Private Room, Sh… Kallang Privat…   110 (31000.55 32666.66)
 9 1-Pax Small Private Room, Sh… Kallang Privat…    78 (31000.55 32666.66)
10 1950s cozy window singlebed … Kallang Privat…    49 (31207.55 32608.05)
# … with 343 more rows, and abbreviated variable names ¹​neighbourhood,
#   ²​room_type
plot(network_sf)

kallang airbnb plot

tmap_mode('plot')
tm_shape(kallang_airbnb_sf) + 
  tm_dots() + 
  tm_shape(network_sf) +
  tm_lines()

all airbnbs plot

tmap_mode('view')
tm_shape(airbnb_sf) + 
  tm_dots() + 
  tm_shape(network_sf) +
  tm_lines()
tmap_mode('plot')

convert to single linestring

network_linestring <- network_sf %>% st_cast("LINESTRING")

save to rds

network_linestring_rds <- write_rds(network_linestring, "data/network_linestring.rds")

preparing lixels objects

lixels <- lixelize_lines(network_linestring, 
                         700, 
                         mindist = 350)

save to rds

lixels_rds <- write_rds(lixels, "data/lixels_rds.rds")

read the rds file

lixels_rds <- read_rds("data/lixels_rds.rds")

generate line centre points

samples <- lines_center(lixels_rds)

write to rds

samples_rds <- write_rds(samples, "data/samples_rds.rds")

read the rds file

samples_rds <- read_rds("data/samples_rds.rds")

perform NetKDE

densities <- nkde(network_linestring, 
                  events = kallang_airbnb_sf,
                  w = rep(1,nrow(kallang_airbnb_sf)),
                  samples = samples_rds,
                  kernel_name = "quartic",
                  bw = 300, 
                  div= "bw", 
                  method = "simple", 
                  digits = 1, 
                  tol = 1,
                  grid_shape = c(1,1), 
                  max_depth = 8,
                  agg = 5, #we aggregate events within a 5m radius (faster calculation)
                  sparse = TRUE,
                  verbose = FALSE)

write to rds

densities_rds <- write_rds(densities, "data/densities_rds.rds")

read the rds file

densities_rds <- read_rds("data/densities_rds.rds")
samples_rds$density <- densities_rds
lixels_rds$density <- densities_rds
# rescaling to help the mapping
samples_rds$density <- samples_rds$density*1000
lixels_rds$density <- lixels_rds$density*1000
tmap_mode('view')
tm_shape(lixels_rds)+
  tm_lines(col="density")+
tm_shape(kallang_airbnb_sf)+
  tm_dots()
kfun_kallang_airbnb <- kfunctions(network_linestring, 
                             kallang_airbnb_sf,
                             start = 0, 
                             end = 1000, 
                             step = 50, 
                             width = 50, 
                             nsim = 50, 
                             resolution = 50,
                             verbose = FALSE, 
                             conf_int = 0.05,
                             agg = 40)
kfun_kallang_airbnb$plotk

for all airbnbs

perform NetKDE for all airbnbs

all_densities <- nkde(network_linestring, 
                  events = airbnb_sf,
                  w = rep(1,nrow(airbnb_sf)),
                  samples = samples_rds,
                  kernel_name = "quartic",
                  bw = 300, 
                  div= "bw", 
                  method = "simple", 
                  digits = 1, 
                  tol = 1,
                  grid_shape = c(1,1), 
                  max_depth = 8,
                  agg = 5, #we aggregate events within a 5m radius (faster calculation)
                  sparse = TRUE,
                  verbose = FALSE)

write to rds

all_densities_rds <- write_rds(all_densities, "data/all_densities_rds.rds")

read rds

all_densities_rds <- read_rds("data/all_densities_rds.rds")
samples_rds$density <- all_densities_rds
lixels_rds$density <- all_densities_rds
# rescaling to help the mapping
samples_rds$density <- samples_rds$density*1000
lixels_rds$density <- lixels_rds$density*1000
tmap_mode('view')
tm_shape(lixels_rds)+
  tm_lines(col="density")+
tm_shape(airbnb_sf)+
  tm_dots()
kfun_all_airbnbs <- kfunctions(network_linestring, 
                             airbnb_sf,
                             start = 0, 
                             end = 1000, 
                             step = 50, 
                             width = 50, 
                             nsim = 50, 
                             resolution = 50,
                             verbose = FALSE, 
                             conf_int = 0.05,
                             agg = 1500)
kfun_all_airbnbs$plotk

network cross K function spNetwork

crossk_kallang <- cross_kfunctions(network_linestring, kallang_airbnb_sf, mrt_sf, start = 0, end = 5000, step = 50, width = 1000, nsim = 50, verbose = FALSE, agg = 100)
crossk_kallang$plotk
crossk_kallang <- cross_kfunctions(network_linestring, mrt_sf, kallang_airbnb_sf, start = 0, end = 5000, step = 50, width = 1000, nsim = 50, verbose = FALSE, agg = 100)
crossk_kallang$plotk

Analysing Marked Point Patterns

make sf containing all 8 different types of point events

attr_mp <- attr_sf
attr_mp$type <- "attraction"
attr_mp <- attr_mp %>% select(c(2,3))

busstop_mp <- busstop_sf
busstop_mp$type <- "busstop"
busstop_mp <- busstop_mp %>% select(c(2,3))

hotel_mp <- hotel_sf
hotel_mp$type <- "hotel"
hotel_mp <- hotel_mp %>% select(c(2,3))

mall_mp <- mall_sf
mall_mp$type <- "mall"
mall_mp <- mall_mp %>% select(c(2,3))

mrt_mp <- mrt_sf
mrt_mp$type <- "mrt"
mrt_mp <- mrt_mp %>% select(c(2,3))

sevenele_mp <- sevenele_sf
sevenele_mp$type <- "seveneleven"
sevenele_mp <- sevenele_mp %>% select(c(2,3))

unis_mp <- unis_sf
unis_mp$type <- "uni"
unis_mp <- unis_mp %>% select(c(2,3))

airbnb_mp <- airbnb_sf
airbnb_mp$type <- "airbnb"
airbnb_mp <- airbnb_mp %>% select(c(5,6))
all_sf <- do.call("rbind", list(attr_mp, busstop_mp, hotel_mp, mall_mp, mrt_mp, sevenele_mp, unis_mp, airbnb_mp))

# all_sf <- do.call("rbind", list(busstop_mp, hotel_mp, mall_mp, mrt_mp, sevenele_mp, unis_mp, airbnb_mp))
all_sf[rowSums(is.na(all_sf))!=0,]
Simple feature collection with 0 features and 1 field
Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
Projected CRS: SVY21 / Singapore TM
[1] type     geometry
<0 rows> (or 0-length row.names)
all_sf[!st_is_empty(all_sf), ]
Simple feature collection with 9979 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 3970.122 ymin: 22858.67 xmax: 48284.56 ymax: 52983.82
Projected CRS: SVY21 / Singapore TM
First 10 features:
         type                  geometry
1  attraction POINT (30007.41 30267.17)
2  attraction POINT (30877.14 31594.08)
3  attraction POINT (29342.34 29382.57)
4  attraction  POINT (29818.66 30598.9)
5  attraction POINT (30043.47 30820.05)
6  attraction POINT (30113.58 30488.32)
7  attraction  POINT (29142.6 29293.89)
8  attraction POINT (30244.32 31323.04)
9  attraction POINT (29165.65 29443.75)
10 attraction POINT (29227.71 29603.72)
all <- as_Spatial(all_sf)

shows that type field is in chr data type, need to change to factor data type

str(all)
Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
  ..@ data       :'data.frame': 9979 obs. of  1 variable:
  .. ..$ type: chr [1:9979] "attraction" "attraction" "attraction" "attraction" ...
  ..@ coords.nrs : num(0) 
  ..@ coords     : num [1:9979, 1:2] 30007 30877 29342 29819 30043 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : NULL
  .. .. ..$ : chr [1:2] "coords.x1" "coords.x2"
  ..@ bbox       : num [1:2, 1:2] 3970 22859 48285 52984
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "coords.x1" "coords.x2"
  .. .. ..$ : chr [1:2] "min" "max"
  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
  .. .. ..@ projargs: chr "+proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +to"| __truncated__
  .. .. ..$ comment: chr "PROJCRS[\"SVY21 / Singapore TM\",\n    BASEGEOGCRS[\"SVY21\",\n        DATUM[\"SVY21\",\n            ELLIPSOID["| __truncated__

convert type field to factor data type

all@data$type <- as.factor(all@data$type)
tmap_mode("plot")
# mpsz_sf
tm_shape(mpsz_sf, bbox=st_bbox(mpsz_sf[mpsz_sf$REGION_N=='CENTRAL REGION',])) +
  tm_borders(alpha = 0.5) +
  tmap_options(check.and.fix = TRUE) +
tm_shape(airbnb_sf) +
  tm_dots(col='red', size = 0.02, alpha = 0.5) +
tm_shape(attr_sf) +
  tm_dots(col='orange', size = 0.02, alpha = 0.5) + 
tm_shape(busstop_sf) +
  tm_dots(col='yellow', size = 0.02, alpha = 0.5) +
tm_shape(sevenele_sf) +
  tm_dots(col='green', size = 0.02, alpha = 0.5) + 
tm_shape(mall_sf) +
  tm_dots(col='blue', size = 0.02, alpha = 0.5) +
tm_shape(hotel_sf) +
  tm_dots(col='purple', size = 0.02, alpha = 0.5) + 
tm_shape(unis_sf) +
  tm_dots(col='black', size = 0.02, alpha = 0.5) +
tm_shape(mrt_sf) +
  tm_dots(col='violet', size = 0.02, alpha = 0.5) +
tm_view(set.zoom.limits = c(10,16)) 

tmap_mode("view")
tm_shape(mpsz_sf) +
  tm_borders(alpha = 0.5) +
  tmap_options(check.and.fix = TRUE) +
tm_shape(all) +
  tm_dots(col = 'type', size = 0.02)
tmap_mode("plot")

converting into ppp

all_ppp <- as(all, "ppp")
plot(all_ppp)
tmap_mode("view")
tm_shape(mpsz_sf) +
  tm_borders(alpha = 0.5) +
  tmap_options(check.and.fix = TRUE) +
tm_shape(all) +
  tm_dots(col = 'type', size = 0.02)
summary(all_ppp)
Marked planar point pattern:  9979 points
Average intensity 7.475023e-06 points per square unit

*Pattern contains duplicated points*

Coordinates are given to 3 decimal places
i.e. rounded to the nearest multiple of 0.001 units

Multitype:
            frequency  proportion    intensity
airbnb           3037 0.304339100 2.274942e-06
attraction        109 0.010922940 8.164922e-08
busstop          5159 0.516985700 3.864480e-06
hotel             422 0.042288810 3.161098e-07
mall              184 0.018438720 1.378299e-07
mrt               562 0.056318270 4.209804e-07
seveneleven       467 0.046798280 3.498182e-07
uni                39 0.003908207 2.921394e-08

Window: rectangle = [3970.12, 48284.56] x [22858.67, 52983.82] units
                    (44310 x 30130 units)
Window area = 1334980000 square units

avoiding duplicated spatial point events using jittering

all_ppp_jit <- rjitter(all_ppp, retry=TRUE, nsim=1, drop=TRUE)

check for duplicates

any(duplicated(all_ppp_jit))
[1] FALSE
summary(all_ppp_jit)
Marked planar point pattern:  9979 points
Average intensity 7.475023e-06 points per square unit

Coordinates are given to 3 decimal places
i.e. rounded to the nearest multiple of 0.001 units

Multitype:
            frequency  proportion    intensity
airbnb           3037 0.304339100 2.274942e-06
attraction        109 0.010922940 8.164922e-08
busstop          5159 0.516985700 3.864480e-06
hotel             422 0.042288810 3.161098e-07
mall              184 0.018438720 1.378299e-07
mrt               562 0.056318270 4.209804e-07
seveneleven       467 0.046798280 3.498182e-07
uni                39 0.003908207 2.921394e-08

Window: rectangle = [3970.12, 48284.56] x [22858.67, 52983.82] units
                    (44310 x 30130 units)
Window area = 1334980000 square units
central = mpsz_sf[mpsz_sf$REGION_N=='CENTRAL REGION',]
kl = mpsz_sf[mpsz_sf$PLN_AREA_N == "KALLANG",]
dc = mpsz_sf[mpsz_sf$PLN_AREA_N == "DOWNTOWN CORE",]
ot = mpsz_sf[mpsz_sf$PLN_AREA_N == "OUTRAM",]
rc = mpsz_sf[mpsz_sf$PLN_AREA_N == "ROCHOR",]

jw = mpsz_sf[mpsz_sf$PLN_AREA_N == "JURONG WEST",]
sb = mpsz_sf[mpsz_sf$PLN_AREA_N == "SEMBAWANG",]
pr = mpsz_sf[mpsz_sf$PLN_AREA_N == "PASIR RIS",]

extracting study area, Kallang, since there is the most airbnbs

kl = mpsz[mpsz@data$PLN_AREA_N == "KALLANG",]
plot(kl, main="Kallang")

converting the spatial point data frame into generic sp format

kl_sp = as(kl, "SpatialPolygons")
str(kl_sp)
Formal class 'SpatialPolygons' [package "sp"] with 4 slots
  ..@ polygons   :List of 9
  .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
  .. .. .. ..@ Polygons :List of 1
  .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 32740 31438
  .. .. .. .. .. .. ..@ area   : num 2434594
  .. .. .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. .. .. ..@ coords : num [1:130, 1:2] 31648 31633 31562 31559 31556 ...
  .. .. .. ..@ plotOrder: int 1
  .. .. .. ..@ labpt    : num [1:2] 32740 31438
  .. .. .. ..@ ID       : chr "52"
  .. .. .. ..@ area     : num 2434594
  .. .. .. ..$ comment: chr "0"
  .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
  .. .. .. ..@ Polygons :List of 1
  .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 32137 32276
  .. .. .. .. .. .. ..@ area   : num 742203
  .. .. .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. .. .. ..@ coords : num [1:63, 1:2] 32246 31916 31744 31699 31682 ...
  .. .. .. ..@ plotOrder: int 1
  .. .. .. ..@ labpt    : num [1:2] 32137 32276
  .. .. .. ..@ ID       : chr "64"
  .. .. .. ..@ area     : num 742203
  .. .. .. ..$ comment: chr "0"
  .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
  .. .. .. ..@ Polygons :List of 1
  .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 30874 32570
  .. .. .. .. .. .. ..@ area   : num 756923
  .. .. .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. .. .. ..@ coords : num [1:103, 1:2] 31422 31387 31371 31357 31342 ...
  .. .. .. ..@ plotOrder: int 1
  .. .. .. ..@ labpt    : num [1:2] 30874 32570
  .. .. .. ..@ ID       : chr "123"
  .. .. .. ..@ area     : num 756923
  .. .. .. ..$ comment: chr "0"
  .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
  .. .. .. ..@ Polygons :List of 1
  .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 29897 32879
  .. .. .. .. .. .. ..@ area   : num 1052033
  .. .. .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. .. .. ..@ coords : num [1:94, 1:2] 30133 30089 29982 29929 29918 ...
  .. .. .. ..@ plotOrder: int 1
  .. .. .. ..@ labpt    : num [1:2] 29897 32879
  .. .. .. ..@ ID       : chr "128"
  .. .. .. ..@ area     : num 1052033
  .. .. .. ..$ comment: chr "0"
  .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
  .. .. .. ..@ Polygons :List of 1
  .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 32338 33019
  .. .. .. .. .. .. ..@ area   : num 410940
  .. .. .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. .. .. ..@ coords : num [1:93, 1:2] 32193 32161 32124 32123 32122 ...
  .. .. .. ..@ plotOrder: int 1
  .. .. .. ..@ labpt    : num [1:2] 32338 33019
  .. .. .. ..@ ID       : chr "129"
  .. .. .. ..@ area     : num 410940
  .. .. .. ..$ comment: chr "0"
  .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
  .. .. .. ..@ Polygons :List of 1
  .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 31579 32929
  .. .. .. .. .. .. ..@ area   : num 839679
  .. .. .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. .. .. ..@ coords : num [1:73, 1:2] 31579 31437 31413 31409 31399 ...
  .. .. .. ..@ plotOrder: int 1
  .. .. .. ..@ labpt    : num [1:2] 31579 32929
  .. .. .. ..@ ID       : chr "130"
  .. .. .. ..@ area     : num 839679
  .. .. .. ..$ comment: chr "0"
  .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
  .. .. .. ..@ Polygons :List of 1
  .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 32110 33766
  .. .. .. .. .. .. ..@ area   : num 735502
  .. .. .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. .. .. ..@ coords : num [1:113, 1:2] 32297 32074 32020 31817 31799 ...
  .. .. .. ..@ plotOrder: int 1
  .. .. .. ..@ labpt    : num [1:2] 32110 33766
  .. .. .. ..@ ID       : chr "147"
  .. .. .. ..@ area     : num 735502
  .. .. .. ..$ comment: chr "0"
  .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
  .. .. .. ..@ Polygons :List of 1
  .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 31176 33828
  .. .. .. .. .. .. ..@ area   : num 1376829
  .. .. .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. .. .. ..@ coords : num [1:123, 1:2] 31231 31277 31306 31307 31338 ...
  .. .. .. ..@ plotOrder: int 1
  .. .. .. ..@ labpt    : num [1:2] 31176 33828
  .. .. .. ..@ ID       : chr "151"
  .. .. .. ..@ area     : num 1376829
  .. .. .. ..$ comment: chr "0"
  .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
  .. .. .. ..@ Polygons :List of 1
  .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 31297 31873
  .. .. .. .. .. .. ..@ area   : num 235998
  .. .. .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. .. .. ..@ coords : num [1:72, 1:2] 31452 31494 31587 31654 31640 ...
  .. .. .. ..@ plotOrder: int 1
  .. .. .. ..@ labpt    : num [1:2] 31297 31873
  .. .. .. ..@ ID       : chr "277"
  .. .. .. ..@ area     : num 235998
  .. .. .. ..$ comment: chr "0"
  ..@ plotOrder  : int [1:9] 1 8 4 6 3 2 7 5 9
  ..@ bbox       : num [1:2, 1:2] 29225 30729 33809 34739
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "x" "y"
  .. .. ..$ : chr [1:2] "min" "max"
  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
  .. .. ..@ projargs: chr "+proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +to"| __truncated__
  .. .. ..$ comment: chr "PROJCRS[\"SVY21 / Singapore TM\",\n    BASEGEOGCRS[\"SVY21\",\n        DATUM[\"SVY21\",\n            ELLIPSOID["| __truncated__
  ..$ comment: chr "TRUE"

creating owin object

kl_owin = as(kl_sp, "owin")
str(kl_owin)
List of 5
 $ type  : chr "polygonal"
 $ xrange: num [1:2] 29225 33809
 $ yrange: num [1:2] 30729 34739
 $ bdry  :List of 9
  ..$ :List of 2
  .. ..$ x: num [1:129] 31648 31661 31773 32097 32186 ...
  .. ..$ y: num [1:129] 30744 30747 30767 30827 30838 ...
  ..$ :List of 2
  .. ..$ x: num [1:62] 32246 32310 32333 32353 32365 ...
  .. ..$ y: num [1:62] 32061 32093 32103 32110 32113 ...
  ..$ :List of 2
  .. ..$ x: num [1:102] 31422 31429 31452 31437 31429 ...
  .. ..$ y: num [1:102] 32135 32141 32159 32193 32223 ...
  ..$ :List of 2
  .. ..$ x: num [1:93] 30133 30156 30175 30200 30214 ...
  .. ..$ y: num [1:93] 32537 32578 32630 32700 32727 ...
  ..$ :List of 2
  .. ..$ x: num [1:92] 32193 32213 32222 32316 32446 ...
  .. ..$ y: num [1:92] 32561 32572 32578 32650 32741 ...
  ..$ :List of 2
  .. ..$ x: num [1:72] 31579 31648 31768 31819 31843 ...
  .. ..$ y: num [1:72] 32380 32396 32426 32438 32443 ...
  ..$ :List of 2
  .. ..$ x: num [1:112] 32297 32320 32338 32346 32352 ...
  .. ..$ y: num [1:112] 33255 33275 33300 33321 33348 ...
  ..$ :List of 2
  .. ..$ x: num [1:122] 31231 31231 31229 31212 31174 ...
  .. ..$ y: num [1:122] 34739 34739 34642 34538 34438 ...
  ..$ :List of 2
  .. ..$ x: num [1:71] 31452 31429 31422 31387 31371 ...
  .. ..$ y: num [1:71] 32159 32141 32135 32097 32081 ...
 $ units :List of 3
  ..$ singular  : chr "unit"
  ..$ plural    : chr "units"
  ..$ multiplier: num 1
  ..- attr(*, "class")= chr "unitname"
 - attr(*, "class")= chr "owin"
kl_owin <- mpsz_sf[mpsz_sf$PLN_AREA_N == "KALLANG",] %>%
  as('Spatial') %>%
  as('SpatialPolygons') %>%
  as('owin')

combining all spatial event points and kallang study area

all_kl_ppp = all_ppp_jit[kl_owin]
summary(all_kl_ppp)
Marked planar point pattern:  548 points
Average intensity 6.383449e-05 points per square unit

Coordinates are given to 3 decimal places
i.e. rounded to the nearest multiple of 0.001 units

Multitype:
            frequency  proportion    intensity
airbnb            344 0.627737200 4.007128e-05
attraction          1 0.001824818 1.164863e-07
busstop           112 0.204379600 1.304646e-05
hotel              47 0.085766420 5.474856e-06
mall                7 0.012773720 8.154041e-07
mrt                21 0.038321170 2.446212e-06
seveneleven        15 0.027372260 1.747294e-06
uni                 1 0.001824818 1.164863e-07

Window: polygonal boundary
9 separate polygons (no holes)
           vertices    area relative.area
polygon 1       129 2434590        0.2840
polygon 2        62  742203        0.0865
polygon 3       102  756923        0.0882
polygon 4        93 1052030        0.1230
polygon 5        92  410940        0.0479
polygon 6        72  839679        0.0978
polygon 7       112  735502        0.0857
polygon 8       122 1376830        0.1600
polygon 9        71  235998        0.0275
enclosing rectangle: [29224.84, 33809.38] x [30728.65, 34738.73] units
                     (4585 x 4010 units)
Window area = 8584700 square units
Fraction of frame area: 0.467
plot(all_kl_ppp,cols=c("red", "orange", "yellow", "brown", "green", "violet", "blue", "black"))

all_Kcross <- Kcross(all_kl_ppp, i="airbnb", j="seveneleven",
                     correction='border')
plot(all_Kcross)
all_Kcross.csr <- envelope(all_kl_ppp, Kcross, i="airbnb", j="seveneleven", correction='border', nsim=999)
Generating 999 simulations of CSR  ...
1, 2, 3, ......10.........20.........30.........40.........50.........60........
.70.........80.........90.........100.........110.........120.........130......
...140.........150.........160.........170.........180.........190.........200....
.....210.........220.........230.........240.........250.........260.........270..
.......280.........290.........300.........310.........320.........330.........340
.........350.........360.........370.........380.........390.........400........
.410.........420.........430.........440.........450.........460.........470......
...480.........490.........500.........510.........520.........530.........540....
.....550.........560.........570.........580.........590.........600.........610..
.......620.........630.........640.........650.........660.........670.........680
.........690.........700.........710.........720.........730.........740........
.750.........760.........770.........780.........790.........800.........810......
...820.........830.........840.........850.........860.........870.........880....
.....890.........900.........910.........920.........930.........940.........950..
.......960.........970.........980.........990........ 999.

Done.
plot(all_Kcross.csr, xlab="distance(m)", xlim=c(0,500))

L cross

all_Lcross <- Lcross(all_kl_ppp, i="airbnb", j="seveneleven", correction='border')
plot(all_kl_ppp, . -r ~ r, 
     xlab = "distance(m)", 
     xlim=c(0, 500))

all_Lcross.csr <- envelope(all_kl_ppp, Lcross, i="airbnb", j="seveneleven", correction='border', nsim=999)
Generating 999 simulations of CSR  ...
1, 2, 3, ......10.........20.........30.........40.........50.........60........
.70.........80.........90.........100.........110.........120.........130......
...140.........150.........160.........170.........180.........190.........200....
.....210.........220.........230.........240.........250.........260.........270..
.......280.........290.........300.........310.........320.........330.........340
.........350.........360.........370.........380.........390.........400........
.410.........420.........430.........440.........450.........460.........470......
...480.........490.........500.........510.........520.........530.........540....
.....550.........560.........570.........580.........590.........600.........610..
.......620.........630.........640.........650.........660.........670.........680
.........690.........700.........710.........720.........730.........740........
.750.........760.........770.........780.........790.........800.........810......
...820.........830.........840.........850.........860.........870.........880....
.....890.........900.........910.........920.........930.........940.........950..
.......960.........970.........980.........990........ 999.

Done.
plot(all_Lcross.csr, . -r ~ r, xlab="distance(m)", xlim=c(0,500))

LCLQ

prepare nb list

nb <- include_self(
  st_knn(st_geometry(all_sf), 6))

weights list

wt <- st_kernel_weights(nb, 
                        all_sf, 
                        "gaussian", 
                        adaptive = TRUE)

vector list

Airbnb <- all_sf %>%
  filter(type == "airbnb")
A <- Airbnb$type
SevenEleven <- all_sf %>%
  filter(type == "seveneleven")
B <- SevenEleven$type

compute LCLQ

LCLQ <- local_colocation(A, B, nb, wt, 49)

joining output table

LCLQ_all <- cbind(all_sf, LCLQ)

plot LCLQ values

tmap_mode("view")
tm_shape(mpsz_sf) +
  tm_polygons() +
tm_shape(LCLQ_all)+ 
  tm_dots(col = "X7.Eleven",
             size = 0.01,
             border.col = "black",
             border.lwd = 0.5) +
  tm_view(set.zoom.limits = c(12, 16))